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In  the  present  study  we  report  on  numerical  investigations  into  the  effects  of  compression  on  the  per¬ 
formance  of  a  unit  cell.  The  focus  of  this  study  is  how  the  transport  properties  of  the  gas  diffusion  layer 
(GDL)  material,  specifically  porosity  and  permeability,  affect  numerical  predictions  of  cell  performance. 
Experimental  data  of  porosity  and  permeability  of  uncompressed  and  compressed  GDLs  were  obtained 
using  a  porometer,  and  used  in  numerical  simulations.  A  3D  model  with  two  parallel  channels  and  an 
membrane  electrode  assembly  (MEA)  is  constructed  for  the  calculations.  Three  different  configurations 
of  transport  properties  were  tested,  i.e.  uniform  uncompressed  GDL  properties,  uniform  compressed  GDL 
properties,  and  non-homogeneous  GDL  properties.  It  is  found  that  the  non-homogeneous  case  shows 
noticeable  differences  in  predicted  cell  performance.  For  the  non-homogenous  case,  simulations  with  a 
pressure  difference  between  two  cathode  channels  were  carried  out  to  gain  insight  into  the  effect  of  cross¬ 
channel  flow  on  the  overall  prediction  of  cell  performance.  We  found  that  the  cross-channel  flow  changes 
local  current  density  distribution  primarily  on  the  high-pressure  channel.  The  present  study  demonstrates 
the  importance  of  the  proper  use  of  transport  properties  for  the  compressed  portion  of  the  GDL. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  gas  diffusion  layer  (GDL),  typically  made  of  carbon  paper 
or  carbon  cloth  and  coated  with  Teflon  and  binder  materials,  is  a 
key  component  in  a  proton  exchange  membrane  fuel  cell  (PEMFC) 
because  of  its  central  role  in  transport  processes  in  the  unit  cell.  The 
GDL  has  multiple  functions:  it  serves  as  a  support  for  the  catalyst 
layer  and  the  membrane,  and  as  a  conductor  for  electricity  and  heat, 
and  provides  pathways  for  species  transport. 

Fig.  1  shows  a  cross-section  of  the  unit  cell  with  the  mem¬ 
brane  electrode  assembly  (MEA)  sandwiched  between  two  flow 
field  plates.  The  GDL  is  in  contact  with  the  bipolar  plate  on  the 
one  side  and  with  the  catalyst  layer  on  the  other.  In  typical  plate- 
and-frame  designs  of  PEMFC,  the  bipolar  plate  has  flow  channels  for 
gas  distribution  and  land  area,  termed  ‘ribs’,  that  remain  in  contact 
with  the  GDL  for  electric  conduction.  In  order  to  achieve  optimum 
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performance,  a  uniform  distribution  of  gas  species  on  the  catalyst 
layer,  and  minimum  electrical  and  thermal  resistance  are  required. 
When  unit  cells  are  assembled  into  a  stack,  a  compression  force 
is  applied  to  the  entire  stack  to  minimize  gas  leakage  through  the 
unit  cell.  However,  these  requirements  often  conflict,  which  calls 
for  careful  design  of  the  bipolar  plate. 

The  characteristics  of  the  GDL  as  a  porous  medium  have  sig¬ 
nificant  impact  on  cell  performance.  Important  factors  include 
volumetric  properties  such  as  porosity,  permeability,  tortuosity, 
thermal  and  electrical  conductivity,  and  surface  properties  such  as 
wettability  and  roughness.  The  effects  of  these  properties  on  cell 
performance  have  been  published  in  literature  [1-12].  In  general, 
mass  transport  of  gaseous  species  driven  by  diffusion  is  affected  by 
porosity  and  tortuosity,  whereas  mass  transport  due  to  pressure  dif¬ 
ference  is  affected  by  permeability,  e.g.  cross-channel  gas  flow  and 
liquid  water  flow  from  catalyst  layer  to  gas  channel  [13,14].  Ther¬ 
mal  and  electrical  conduction  are  comparable  because  of  the  fact 
that  both  transport  processes  take  place  in  similar  (solid)  phase, 
although  for  thermal  conduction  the  binder  and  gas  in  the  pores 
also  make  a  contribution,  in  addition  to  that  made  via  the  carbon 
fibers.  The  thermal  and  electrical  contact  resistances  of  the  GDL 
with  the  bipolar  plate  are  sensitive  to  the  surface  roughness  of  the 
GDL.  Transport  of  liquid  water  in  the  GDL  has  been  the  focus  of 
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Fig.  1.  Cross-section  area  of  single  fuel  cell. 


numerous  papers,  e.g.  Refs.  [15-17].  The  transport  of  liquid  water 
is  mostly  dependent  on  the  wettability  of  the  GDL  materials  and 
the  connectivity  of  the  pore  network.  Since  both  liquid  water  and 
gaseous  species  occupy  the  same  pore  volume,  transport  of  liquid 
water  has  a  profound  impact  on  the  mass  transfer  of  gas  within  the 
GDL.  The  mass,  heat  and  electron  transport  processes  in  the  GDL  are 
coupled  mainly  due  to  the  rib/channel  configuration  that  causes  a 
skewed  distribution  of  mass,  heat  and  electron  flows  [18,19].  The 
problem  is  further  complicated  when  deformation  occurs  in  the 
GDL  due  to  compression  force  exerted  on  the  stack.  The  material 
properties  of  both  the  rib  area  and  the  channel  area  of  the  GDL  may 
change  dramatically  due  to  such  deformation. 

Models  and  simulation  tools  have  been  developed  to  facilitate 
analysis  of  the  complex  problems  involved  in  understanding  the 
GDL,  e.g.  Ref.  [20].  Uniform  material  properties  for  the  GDL  are 
adopted  by  most  authors  [21-26].  In  these  numerical  simulations, 
constant  and  uniform  transport  properties  were  set  for  the  GDL 
material  under  the  ribs  and  the  channels.  Such  cases  are  referred  to 
as  homogeneous  cases  in  the  present  study.  Apparently,  the  homo¬ 
geneous  case  cannot  predict  the  effects  of  compression  force  on  cell 
performance.  There  have  been  few  systematic  investigations  of  the 
effects  of  non-uniform  transport  properties  of  the  GDL  due  to  com¬ 
pression.  Gurau  et  al.  [27]  developed  a  one-dimensional  half-cell 
model  that  considered  non-uniform  GDL  properties  to  account  for 
the  fact  that  the  pores  might  be  partially  filled  with  liquid  water. 
Chu  et  al.  [28]  developed  a  one-dimensional  half-cell  model  and 
studied  the  impact  of  non-uniformity  GDL  porosity  with  four  differ¬ 
ent  continuous  functions  of  the  position  (constant,  linear,  convex, 
and  concave  exponential  function).  Roshandel  et  al.  [29]  reported 
a  two-dimensional  model  that  simulates  the  porosity  change  of 
the  GDL  after  compression.  They  used  a  composition  function  in 
the  form  of  Sin2n(x)  to  describe  the  effect  on  cell  performance  of 
a  change  in  the  porosity  of  the  GDL.  Zhou  et  al.  [30]  investigated 
the  effect  of  clamping  force  on  the  performance  of  PEMFC  with 
an  interdigitated  gas  distributor,  considering  the  interfacial  contact 
resistance,  the  non-uniform  porosity  distribution  of  the  GDL,  and 
the  GDL  deformation.  In  a  recent  work,  Zhou  and  Wu  [31  ]  included 
transport  of  liquid  water  in  the  model  and  reported  simulation 
results  in  a  2D  configuration.  The  effects  of  non-uniform  compres¬ 
sion  of  the  GDL  under  the  channel/rib  structure  of  flow-filled  plate 
on  the  temperature  distribution  in  the  PEMFC  is  studied  in  Hotti- 
nen  and  Himanen  [32].  Recent  work  of  Nitta  et  al.  [33]  and  Hottinen 
et  al.  [34]  reported  more  thorough  experimental  and  numerical 
investigations  of  this  issue. 

In  the  present  study  we  report  on  the  results  of  numerical 
investigations  when  the  change  of  the  transport  properties  of 
GDL  due  to  compression  force  is  considered.  Different  porosity 
and  permeability  measured  from  compressed  and  uncompressed 


GDL  material  are  assigned  to  the  rib  and  channel  areas  of  the 
GDL,  respectively.  This  case  is  termed  the  non-homogeneous  case 
in  the  following  discussion.  The  objective  of  the  present  study 
is  to  investigate  the  effects  of  compression  force  on  cell  perfor¬ 
mance.  A  three-dimensional  computational  model  using  properties 
and  parameters  validated  with  experimental  data  is  constructed 
to  assess  the  sensitivity  of  predicted  results  when  realistic,  non- 
uniform  GDL  properties  are  used. 

2.  Characterization  of  the  GDL 

Fig.  1  shows  a  cross-section  area  of  a  unit  cell  sandwiched  by 
two  pieces  of  plates  as  the  compression  fixture.  The  GDL  under 


Fig.  2.  SEM  photo  of  the  compressed  and  uncompressed  GDL  as  (a)  obersvered  zone; 
(b)  500  |jim;  (c)  10|jim. 
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Fig.  3.  Pore  size  distribution  (left:  uncompressed  GDL;  right:  compressed  GDL). 


the  rib  area  appears  to  be  thinner  than  the  portion  in  the  channel 
area.  Since  the  GDL  serves  multiple  purposes  of  transport  (reactant 
gases,  product  water,  electrons,  and  heat),  characterization  of  the 
GDL  material  in  both  regions  will  help  modeling  and  simulations 
of  the  transport  in  the  GDL.  In  the  present  study  we  measure  the 
permeability  and  porosity  of  the  compressed  (under  the  rib)  and 
uncompressed  (under  the  channel)  zone  of  the  GDL  and  use  the 
measured  properties  in  numerical  simulations  for  prediction  of  cell 
performance  due  to  compression. 

The  bipolar  plate  and  gasket  on  either  side  of  the  MEA  are 
bolted  together  under  significant  clamping  force.  This  enables  opti¬ 


mum  gastight  during  fuel  cell  assembly.  GDL  deformation  depends 
on  the  thickness,  form,  and  material  of  the  gasket  and  clamp¬ 
ing  force  during  fuel  cell  assembly.  However,  compression  force 
is  known  to  change  the  porous  medium  of  the  electrode  and  GDL. 
The  magnitude  of  the  compression  force  affects  directly  diffusion 
and  permeability  of  the  reactant  gases  in  the  compressed  GDL.  A 
high  compression  force  increases  mass  transfer  resistance  in  the 
compressed  GDL,  which  has  lower  porosity  than  uncompressed 
GDL.  Conversely,  a  low  compression  force  increases  the  contact 
resistance  between  the  GDL  and  the  bipolar  plate,  as  well  as  gas 
leakage. 
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Fig.  4.  Computational  domain  (a)  and  mesh  (b)  in  theXY  plane. 


A  scanning  electron  microscope  (SEM,  JSM-5610LV;  JEOL)  was 
employed  to  observe  the  changes  in  the  GDL  surface  microstruc¬ 
ture  compressed  by  the  rib  area  of  the  bipolar  plate.  The  procedure 
of  applying  compression  force  to  a  fuel  cell  assembly  is  illustrated 
in  Ref.  [35].  Fig.  2(a)  shows  a  cross-section  of  a  single  cell,  showing 
two  gas  channels  and  one  rib  in  the  image.  A  close-up  of  the  GDL 
surface  after  compression  by  the  rib  is  shown  in  Fig.  2(b)  and  (c), 
which  depict  a  top  view  of  the  GDL  surface  for  the  region  marked 
with  a  dotted  box  in  Fig.  2(a).  It  is  clear  that,  in  the  zone  under 
compression  by  the  rib,  there  are  two  cracks  on  the  GDL  surface. 
After  the  fuel  cell  is  assembled,  the  physical  characteristics  of  the 
compressed  area  of  the  GDL  become  different  than  those  of  the 
uncompressed  area.  The  GDL  surface  of  the  area  that  is  exposed  to 
the  rib  clearly  contains  more  cracks.  In  Fig.  2(a),  the  GDL  exposed 
to  the  gas  channel  is  not  under  compression  and  only  a  small  por¬ 
tion  protrudes  into  the  channel.  Under  the  rib  of  the  bipolar  plate, 
the  GDL  is  thinner  than  the  channel  portion.  Fig.  2(b)  shows  the 
GDL  with  microporous  coating  on  one  side.  Two  crack  lines  may 


be  seen  on  the  microporous  surface.  Fig.  2(c)  shows  a  close-up  of 
the  area  near  a  crack  line  with  a  depth  of  approximately  10-20  [xm. 
In  the  present  study,  carbon  black  with  30wt.%  PTFE  (Teflon  30J, 
Dupont™)  was  coated  on  the  carbon  cloth  (CPW-003  Textron)  to 
make  the  GDL.  The  properties  of  uncompressed  GDL  were  measured 
and  used  to  represent  the  uncompressed  zone  in  an  assembled  fuel 
cell,  i.e.  the  portion  of  GDL  facing  the  gas  channel.  For  the  properties 
representing  the  material  under  the  rib,  data  was  collected  using 
the  portion  of  GDL  material  compressed  by  the  solid  rib  area  during 
fuel  cell  assembly. 

The  permeability  data  of  the  compressed  and  uncompressed 
GDL  were  measured  by  a  capillary  flow  porometer  (Porous  Mate¬ 
rial,  Inc.  CFP-1100-AEX).  The  porosity  of  the  compressed  and  the 
uncompressed  GDL  was  measured  by  a  water  intrusion  porosime- 
ter  (Porous  Material,  Inc.),  which  measured  the  intrusion  volume 
of  water  into  the  hydrophobic  pores  of  the  material  as  a  function  of 
pressure  applied  on  water.  Information  of  pore  size  distribution  was 
deduced  from  the  pressure  data.  In  the  present  study,  an  average 
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Fig.  5.  Predicted  current  density  profiles  vs.  number  of  grid. 


of  three  to  five  specimens  were  measured  and  data  were  calcu¬ 
lated  to  obtain  a  mean  average  porosity  and  permeability  value. 
Measured  data  with  these  samples  were  within  5%  error.  The  prin¬ 
ciples  for  these  measurement  techniques  and  data  reduction  can 
be  found  in  Refs.  [36,37].  The  sampling  principle  involved  compar¬ 
ing  the  characteristics  of  the  dry  and  wet  samples.  A  wetting  liquid 
is  first  allowed  to  spontaneously  fill  the  pores  in  the  sample  and 
a  non-reacting  gas  is  then  introduced  to  displace  the  liquid  from 
the  pores.  The  gas  pressure  and  flow  rates  through  the  dry  and 
wet  samples  are  accurately  recorded.  The  gas  pressure  required  for 
removing  liquid  from  the  pores  and  causing  gas  to  flow  is  given  by 

(D 

where  D  is  the  pore  diameter,  y  is  the  surface  tension  of  liquid,  0  is 
the  contact  angle  of  liquid,  and  Pis  the  pressure  difference  between 
liquid  and  gas.  A  pore  size  distribution  function, /is  defined  as 

,  d[100  x  (Fw/Fd)] 

}  d  D  1  j 

Table  1 

Flow  channel  dimensions  and  inlet  conditions  used  in  simulations 


Parameter  Value 


Channel  length  (mm)  50 

Channel  width  (mm)  1 .5 

Channel  depth  (mm)  1 

Solid  rib  area  width  (mm)  2 

Gas  diffusion  layer  thickness  (mm)  0.3 

Membrane  thickness  (mm)  0.05 

Catalyst  layer  thickness  (mm)  0.03 

Membrane  porosity  0.28 

Membrane  ionic  conductivity  (S  itt1  )  17 

Electrode  electronic  conductivity  (S  m-1 )  570 

Membrane  thermal  conductivity  (W(mK)-1 )  0.455 

Anode  mass  flow  rate  (kgs-1)  1.79  x  10-7 

Cathode  mass  flow  rate  (kgs-1)  1.28  x  10-6 

Hydrogen  stoichiometric  ratio  1.2 

Air  stoichiometric  ratio  2 

Anode  hydrogen  mass  fraction  0.174 

Anode  water  vapor  mass  fraction  0.826 

Cathode  oxygen  mass  fraction  0.195 

Cathode  nitrogen  mass  fraction  0.645 

Cathode  water  vapor  mass  fraction  0.16 


The  pore  size  distribution  function  calculated  from  the  variation 
of  the  flow  rate  with  the  differential  pressure  for  through-plane 
flow,  shown  in  Fig.  3.  The  area  under  the  distribution  curve  in  a 
given  pore  diameter  range  yields  the  percentage  of  flow  through 
pores  in  that  size  range.  In  the  present  study,  the  porosity  mea¬ 
sured  was  0.42  and  0.64  for  the  compressed  and  uncompressed  GDL, 
respectively.  The  permeability  was  measured  to  6.42  x  10-13  and 
3.64  x  10-11  m2  for  the  compressed  and  uncompressed  GDL,  respec¬ 
tively.  These  values  were  used  in  the  computer  simulations  to  be 
discussed  later. 

3.  Mathematical  formulation  and  numerical  simulations 

3.1.  Assumptions 

The  reaction  of  the  fuel  cell  is  complex,  including  electrochem¬ 
ical  reactions,  hydrodynamics,  phase  changes,  heat  transfer,  and 
mass  transfer.  In  order  to  carry  out  the  numerical  simulation  and 
analysis,  effective  simplifications  and  assumptions  must  be  made. 
The  following  are  the  key  assumptions  of  this  model: 

(1 )  The  gas  mixture  is  a  perfect  gas. 

(2)  Steady  state. 

(3)  Constant  cell  temperature. 

(4)  Laminar  flow  and  incompressibility,  because  the  gradient  is  of 
very  small  pressure. 

(5)  The  porous  medium  is  isotropic  and  homogeneous. 

(6)  Single  phase  flow  that  neglects  the  existence  of  the  liquid  water. 

(7)  Deformation  of  the  porous  medium  is  ignored  (shrinks  and 
expands). 

(8)  The  reactant  is  not  permeable  through  the  proton  exchanges 
the  membrane. 

(9)  The  Butler-Volmer  equation  controls  electro-chemical  dynam¬ 
ics. 

A  commercial  software,  CFD-ACE+,  is  used  in  the  present  study. 
Mathematical  formulation  of  the  fuel  cell  module  of  this  software 
can  be  found  in  Mazumder  and  Cole  [38]. 

3.2.  Computational  domain  and  mesh 

A  computational  domain  of  a  double-channel  configuration 
is  adopted  to  simulate  transport  in  a  single  cell  and  to  inves¬ 
tigate  the  effect  of  compression  pressure  on  cell  performance. 
Fig.  4  shows  the  geometry  and  computational  domain  of  the 
double-channel.  The  dimensions  of  the  flow  channel  are  50  mm 
(L)  x  1.5  mm  (W)  x  1  mm  (FI).  The  dimensions  of  the  GDL  are  50  mm 
(L)  x  3.5  mm  (W)  x  0.3  mm  (H). 

Verification  of  grid  independence  for  the  numerical  solutions 
was  performed  to  confirm  the  accuracy  of  the  computational 
results.  Fig.  5  shows  the  relationship  between  the  grid  number 
and  current  density  along  with  the  channel  (Z-axis).  Local  current 
density  along  the  channel  and  the  average  plane  current  den¬ 


iable  2 

The  GDL  characteristics  used  in  the  simulations 


Properties 

Case  type 

Homogeneous(c) 

Homogeneous(u) 

Non-homogeneous 

Rib  area 

Porosity 

Permeability  (m2) 

0.42 

6.42  x  10-13 

0.64 

3.64  x  10-11 

0.42 

6.42  x  10-13 

Channel  area 

Porosity 

Permeability  (m2) 

0.42 

6.42  x  10-13 

0.64 

3.64  x  10-11 

0.64 

3.64  x  10-11 
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Fig.  6.  Predicted  polarization  curve  and  power  for  homogeneous(c),  homoge¬ 
neous^),  and  non-homogeneous  cases. 

sity  were  computed  to  evaluate  the  grid  independence:  when  the 
grid  number  reaches  242,550  grids,  the  variation  of  the  current 
density  trends  towards  stability  at  Z= 0.0125  m.  The  variation  of 
the  current  density  is  not  influenced  by  further  increases  in  the 


Anode 


Mass  Fraction  of  I  I  I  I  I  I  I  I  1  I  I 

Hydrogen  0.004  0.020  0.036  0.052  0.068  0.084  O.lOO  0.116  0.132  0.148  0.164 


Cathode 


0.05 


Mass  Fraction  of 
Oxygen 


. 


0.000  0.020  0.040  0.060  0.080  0.100  0.120  0.140  0.160  0.180  0.200 


grid  number.  The  results  of  the  other  positions  are  similar  too. 
There  is  0.42%  of  difference  in  predicted  results  between  the  cases 
when  242,550  and  338,550  cells  were  used  in  the  domain,  respec¬ 
tively.  The  computational  domain  with  242,550  cells  was  used 
throughout  the  simulations  reported  in  this  paper  unless  otherwise 
noted. 


3.3.  Boundary  conditions  and  properties 

Table  1  lists  the  gas  composition  and  channel  dimensions  used 
in  the  computations.  The  anode  gas  mixture  contains  hydrogen  and 
water  vapor  and  the  cathode  gas  mixture  contains  oxygen,  nitro¬ 
gen,  and  water  vapor.  Mass  flow  rate  with  known  gas  composition 
is  prescribed  for  all  channel  inlets.  Fixed  pressures  are  set  for  the 
channel  outlets.  Table  2  lists  the  GDL  properties  used  in  the  sim¬ 
ulations.  Two  homogeneous  cases  are  studied,  the  homogeneous(c) 
case,  which  uses  properties  of  a  compressed  GDL  material  in  the 
calculation,  and  the  homogeneous(u)  case,  which  uses  properties 
of  a  raw  GDL  material.  For  the  non-homogenous  case,  porosity  and 
permeability  data  were  used  in  their  corresponding  zone  in  the 
computational  domain.  For  the  non-homogeneous  case,  two  outlet 
pressure  conditions  are  tested  to  investigate  the  effect  due  to  cross¬ 
channel  flow.  The  baseline  case  has  identical  pressure  for  all  outlets, 
i.e.  PA1  =  PA2  =  Pci  =  Pc 2  =  100  kPa.  In  the  different  pressure  case,  only 
the  cathode  outlet  pressure  is  varied  from  PCi  =  101  to  103  kPa 
while  all  other  outlets  are  kept  the  same  pressure  as  the  baseline 
case. 


Anode 


Mass  Fraction  of 

Hydrogen  0.004  0.020  0.036  0.052  0.068  0.084  0.100  0.116  0.132  0.148  0.164 


Oxygen  o.ooo  0.020  0.040  0.060  o.oso  0.100  0.120  o.uo  0.I60  o.iso  0.200 


Fig.  7.  Predicted  mass  fraction  of  reactant  gases  for  two  different  GDL  cases:  (a)  homogeneous(c);  (b)  non-homogeneous. 
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Channel  width  (in) 


Fig.  8.  Predicted  velocity  vector  and  mass  fraction  of  water  vapor  (Z=  0.025  mm) 
left-half:  homogenous(c);  right-half:  non-homogeneous. 


4.  Results  and  discussion 

The  variation  of  porosity  and  permeability  in  the  GDL  is  expected 
to  primarily  affect  the  transport  of  reactant  gases  in  the  MEA.  When 
the  pressure  of  the  two  adjacent  gas  channels  is  the  same,  the  trans¬ 
port  of  gas  species  is  primarily  dominated  by  mass  diffusion,  which 
is  related  to  GDL  porosity,  because  convective  flow  due  to  pres¬ 
sure  difference  between  the  channels  is  minimal.  The  effect  due 
to  variable  permeability  becomes  important  when  there  exists  a 
pressure  differential  between  adjacent  gas  channels.  In  Section  4.1 , 
numerical  investigations  on  the  effect  of  GDL  properties  on  pre¬ 
dicted  cell  performance  carried  out  for  gas  channels  with  the  same 
inlet  pressure  are  discussed.  In  Section  4.2,  numerical  results  for  the 
non-homogeneous  case  with  different  inlet  pressure  are  compared. 

4 A.  Homogenous  versus  non-homogeneous  GDL  properties 

Fig.  6  compares  the  predicted  polarization  curves  for  the  homo- 
geneous(c),  homogeneous(u),  and  non-homogeneous  cases.  The  inlet 
pressure  of  adjacent  gas  channels  are  set  the  same.  For  cell  voltage 
above  0.9  V,  one  can  see  that  all  polarization  curves  are  very  close 
because  the  overpotential  in  this  regime  is  dictated  by  the  activation 
loss,  therefore  ohmic  loss  and  potential  loss  due  to  mass  transport 
are  negligible.  For  cell  voltage  below  0.4  V,  the  polarization  curves 
show  effects  due  to  mass  transport  limitation,  which  eventually 
results  in  an  abrupt  drop  of  cell  potential  at  higher  current  den¬ 
sity.  The  limiting  current  density  extrapolated  from  these  curves  for 
the  homogeneous(c),  homogeneous(u),  and  non-homogeneous  cases 
is  0.995,  1.20,  and  1.175  A  cm-2,  respectively.  It  is  noted  that  for 


both  the  homogeneous(u)  and  non-homogeneous  cases,  the  porosity 
immediately  under  the  gas  channel  region  is  the  same.  Comparing 
the  limiting  current  density  value  for  these  cases,  one  can  see  that 
apparently  the  porosity  for  the  portion  of  the  GDL  underneath  the 
gas  channel  has  more  impact  on  the  mass  transfer. 

Fig.  7(a)  and  (b)  compare  the  distribution  of  mass  fraction  of 
reactant  gases  for  homogeneous(c)  and  non-homogeneous  cases  at 
0.4  V.  Computational  results  of  zone  A  of  Fig.  4(b)  are  shown  in  these 
figures.  The  GDL  region  underneath  the  channel  area  has  a  higher 
porosity  for  the  non-homogeneous  case,  therefore  the  reactant  gas 
sees  less  resistance.  The  results  show  an  overall  higher  consump¬ 
tion  rate  of  the  reactant  gases,  hence  a  higher  current  density,  along 
flow  channel  for  the  non-homogeneous  case. 

The  velocity  vectors  on  the  X-Y  plane  and  contours  of  water 
mass  fraction  at  Z=  0.025  m  for  the  homogeneous(c)  and  non- 
homogeneous  cases  are  shown  in  Fig.  8.  Computational  results  of 
zone  B  of  Fig.  4(b)  are  shown  in  this  figure.  Because  for  both  cases 
the  computational  results  are  symmetrical  with  respect  to  the  rib 
centerline,  the  results  are  placed  side  by  side  to  facilitate  compari¬ 
son.  One  can  see  for  the  cathode  side,  the  net  flow  of  the  gas  mixture 
is  from  the  GDL  to  the  gas  channel,  whereas  for  the  anode  side  the 
net  flow  is  from  the  gas  channel  to  the  GDL.  It  is  noted  that  the 
velocity  for  each  species  is  a  result  of  the  fluid  velocity  and  diffu¬ 
sion  velocity,  which  is  expressed  as  a  gradient  of  the  mass  fraction. 
For  both  cases  one  can  see  that  in  the  channel  region  the  velocity 
component  in  the  XT  plane  is  small  because  in  the  gas  flow  the  axial 
direction  (Z-coordinate)  is  dominating.  For  the  flow  under  the  rib  in 
the  GDL,  the  non-homogeneous  case  shows  slightly  lower  velocity 
of  the  flow  in  the  GDL  between  the  portion  under  the  channel  and 
under  the  rib  because  of  their  dissimilar  properties.  More  stratifi¬ 
cation  of  water  concentration  near  this  area  is  also  observed  for  the 
non-homogeneous  case. 

Fig.  9  compares  predicted  current  density  profiles  at  channel 
centerline  (x  =  0.00175)  along  the  flow  channel  for  the  homoge- 
neous(c)  and  non-homogeneous  cases.  The  current  density  shown 
in  Fig.  9  is  averaged  over  the  rib-land  length,  thus  it  is  differ¬ 
ent  from  that  shown  in  Fig.  6,  which  is  averaged  over  the  MEA 
area  of  the  entire  domain.  The  average  current  density  for  the 
homogeneous(c)  and  non-homogeneous  case,  cf.  Fig.  6,  is  993.5  and 
1161.3  mAcrn-2,  respectively,  a  difference  of  14.4%.  A  considerable 
difference  for  Z=  0.025-0.05  m  is  noted.  For  both  cases,  only  the 


Channel  Coordinate  Z(m) 


Fig.  9.  Predicted  current  density  profiles  for  different  GDL  properties. 
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Rib  Channel  Rib 


-  Non-homogeneous 
G —  Homogeneous  (c) 


Fig.  10.  Predicted  current  density  in  Y-direction  for  homogenous(c)  and  non-homogenous  cases. 


portion  of  the  GDL  underneath  the  rib  has  different  porosities, 
i.e.  the  non-homogeneous  case  the  porosity  is  52%  higher  of  that 
for  the  homogeneous(c)  case.  Therefore  the  reactant  gases  diffus¬ 
ing  through  the  GDL  are  greater  for  the  non-homogeneous  case, 
resulting  in  a  higher  current  density.  At  Z=  0.0125  m,  the  predicted 
current  density  of  the  homogeneous(c)  and  non-homogeneous  cases 
are  1216  and  1562  mA  cm-2,  respectively.  AtZ=  0.025  m,  the  current 
density  of  the  homogeneous(c)  and  non-homogeneous  case  are  925 
and  1136  mAcm-2,  respectively.  Both  predictions  demonstrate  that 
higher  porosity  of  the  GDL  results  in  a  better  transport  of  the  reac¬ 
tant  gases  and  consequently  a  higher  current  density  when  mass 
transport  becomes  limited.  At  Z=  0.045  m,  variation  of  the  current 
density  along  the  flow  channel  is  similar  in  both  cases. 

The  distribution  of  current  density  along  cross-section  area  of 
the  flow  channel  calculated  using  the  homogenous(c)  and  non- 
homogenous  cases  is  illustrated  in  Fig.  10  for  three  axial  (Z)  locations. 
Computational  results  of  zone  A  of  Fig.  4(b)  are  shown  in  this  figure. 
For  the  non-homogeneous  case,  the  distribution  of  current  density 
is  a  concave  curve,  shown  in  Fig.  10(a).  With  this  case  the  predicted 
current  density  near  the  border  between  the  rib  and  the  flow  chan¬ 
nel  (3371  mA  cm-2 )  is  higher  than  the  current  density  of  the  central 


position  (2341  mAcnrr2),  cf.  Fig.  10(a),  a  difference  of  about  30.6%. 
(The  maximum  current  density  occurs  near  the  corner  of  the  rib. 
This  is  a  balance  of  mass  transfer  resistance  in  the  GDL  and  the 
electrical  resistance  in  the  GDL  material.)  An  implication  of  the 
highly  non-uniform  current  density  near  the  corner  of  the  rib  is 
that  hot  spots  may  occur  and  damage  the  MEA.  The  current  density 
of  other  position  drops  significantly  to  about  211  mAcm-2  near  the 
symmetry  line  of  the  rib.  Further  downstream  the  flow  channel, 
cf.  Fig.  10(b),  the  profile  of  current  density  resembles  a  bell-shape 
curve.  The  difference  in  the  maximum  current  density  predicted 
with  the  homogeneous(c)  and  the  non-homogeneous  case  is  about 
27.4%  for  at  Z=  0.025  m.  Interestingly,  the  current  profiles  become 
similar  near  the  channel  outlet,  cf.  Fig.  10(c).  In  this  case,  the  gas 
concentration  predicted  by  both  cases  is  similar. 

Fig.  11  compares  predictions  of  oxygen  mass  fraction  at  the 
cathode  catalyst  layer  as  function  of  current  density  for  the  homoge- 
neous(c)  and  non-homogeneous  case.  There  is  a  noticeable  difference 
at  high  current  density  conditions.  Under  high  loads,  the  current 
density  distribution  is  influenced  by  variation  of  oxygen  concentra¬ 
tion  along  the  flow  channel.  When  the  oxygen  mass  fraction  is  0.1, 
the  current  density  of  the  non-homogeneous  case  (438  mAcm-2)  is 
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Fig.  11.  Oxygen  mass  fraction  at  the  cathode  catalyst  layer  as  function  of  current 
density  for  homogeneous(c)  and  non-homogeneous  cases. 

higher  than  that  of  the  homogeneous(c)  case  (375  mA  cm-2 ),  a  16.8% 
difference.  For  oxygen  mass  fraction  at  0.025,  the  predicted  current 
density  for  the  homogenous(c)  and  the  non-homogeneous  case  is  849 
and  1040  mA  cm-2,  respectively,  or  a  22.5%  difference.  The  differ¬ 
ence  in  current  density  is  a  result  of  differing  oxygen  diffusion  rates 
due  to  the  different  GDL  porosity  used  in  the  case. 


Fig.  12.  Predicted  polarization  curve  and  power  of  the  non-homogeneous  case  for 
different  outlet  pressure. 

pressure  channel  (right  channel  in  Fig.  14(b)  is  shifted  towards  the 
rib/channel  boundary.  Flowever,  the  current  density  contours  in  the 
left  channel  remain  fairly  symmetrical,  showing  only  slight  dis¬ 
placement  to  the  left  side.  Notably,  the  current  density  peaks  at 
similar  axial  location  but  the  maximum  value  exceeds  the  baseline 
case,  which  has  no  pressure  difference  in  the  channels. 


4.2.  Cross-channel  flow  due  to  pressure  difference 

Predicted  polarization  curves  for  different  cathode  outlet  pres¬ 
sure  ranging  from  1  to  3kPa  are  shown  in  Fig.  12  with  intervals 
of  1  kPa.  The  outlet  pressure  of  the  anode  channels  is  kept  the 
same.  The  polarization  curves  show  that  cell  performance  increases 
with  increasing  of  cathode  outlet  pressure  difference,  e.g.  at  0.6  V 
the  current  density  for  APc  =  l,  2,  and  31<Pa  is  1084,  1114,  and 
1148  mAcm-2,  respectively.  At  higher  load  conditions,  the  current 
density  increases  slightly  with  increasing  both  cathode  outlet  pres¬ 
sure.  The  reason  for  is  because  of  the  convection  induced  by  the 
pressure  difference  between  two  adjacent  channels  that  introduces 
more  reactant  gas  into  the  rib  area.  Fig.  13  shows  the  velocity 
vector  on  the  X-Y  plane  and  mass  fraction  of  water  vapor  con¬ 
tour  at  Z= 0.025  m  for  the  non-homogeneous  cases  with  a  pressure 
differential  between  two  channels.  Computational  results  of  zone 
B  of  Fig.  4(b)  are  shown  in  this  figure.  Compared  with  the  non- 
homogeneous  case  with  no  pressure  difference  (right-half  portion 
in  Fig.  8),  the  case  with  a  APc  =  2kPa  clearly  shows  a  gas  flow 
through  the  rib  area.  This  cross-channel  flow  helps  to  increase  the 
reactant  concentration  underneath  the  rib  area,  hence  an  increase 
of  limiting  current. 

Fig.  14  compares  the  prediction  of  current  density  distribution 
on  the  cathode  catalyst  layer  for  the  non-homogeneous  case  with¬ 
out  (Fig.  12a)  and  with  (Fig.  12b)  a  pressure  difference  between 
channels.  When  there  is  no  cross-channel  flow,  the  current  den¬ 
sity  distribution  appears  to  be  symmetrical  with  respect  to  the 
centerline  of  the  cathode  channel,  cf.  Fig.  14(a).  The  density  cur¬ 
rent  contours  appear  parabolic  downstream  in  this  case  peaks  at 
X= 0.0165  m  and  from  that  point  on  the  current  density  decays 
towards  the  outlet,  primarily  due  to  depletion  of  oxygen  down¬ 
stream.  When  a  convective  cross-channel  flow  occurs  due  to 
pressure  difference,  the  current  density  distribution  on  the  high- 


Fig.  13.  Predicted  velocity  vector  and  mass  fraction  of  water  vapor  (Z=  0.025  mm) 
in  the  non-homogenous  case  with  A Pc  =  2  kPa. 
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2nd  Channel  area  Rib  area  1st  Channel  area 


2nd  Channel  area  Rib  area  1st  Channel  area 


Fig.  14.  Current  density  distribution  at  the  cathode  catalyst  layer  for  two  different 
non-homogeneous  cases:  (a)  A Pc  =  0  kPa;  (b)  A Pc  =  2  kPa. 

5.  Conclusions 

The  present  study  reports  on  numerical  investigations  on  the 
effect  of  compression  force  on  the  performance  of  PEM  fuel  cells. 
Experimental  data  of  porosity  and  permeability  for  uncompressed 
and  compressed  GDL  are  used  in  the  computation.  The  analysis 
focuses  on  the  transport  phenomenon  of  flow  and  reactant  gases  in 
GDL  regions  underneath  the  rib  and  the  channel  area.  A  3D  numer¬ 
ical  model  consisting  of  two  parallel,  straight  channels  along  with 
an  MEA  and  bipolar  plates  are  constructed  to  study  the  problem. 
A  series  of  numerical  simulations  are  carried  out  to  investigate 
the  distribution  of  GDL  properties  on  cell  performance  and  to 


study  the  effects  of  cross-channel  flow  induced  by  a  pressure  gra¬ 
dient  between  the  gas  channels.  The  following  conclusions  were 
obtained: 

(1)  The  polarization  for  the  cases  tested  {homogeneous(c),  homo¬ 
geneous (u),  and  non-homogeneous  cases)  shows  that  cell 
performance  for  the  configurations  tested  is  controlled  by  the 
transport  under  the  channel  area.  Therefore  the  limiting  current 
density  depends  on  the  porosity  assigned  for  this  portion  of  the 
GDL.  The  predicted  cell  performance  is  homogeneous(u)>  non- 
homogeneous  >  homogeneous(c). 

(2)  In  the  non-homogeneous  case,  a  highly  non-uniform  current 
density  exists  near  the  corner  of  the  rib  where  hot  spots  may 
develop  during  operation.  Therefore,  careful  and  appropriate 
application  of  the  compression  force  as  well  as  rib/channel 
dimensions  and  geometry  are  required  in  order  to  avoid  the 
formation  of  hot  spots  inside  the  MEA. 

(3)  When  the  pressure  difference  between  two  adjacent  channels  is 
high  enough,  a  convective  cross-channel  flow  occurs.  A  predic¬ 
tion  of  the  effect  of  the  cross-channel  flow  shows  that  the  flow 
can  bring  more  reactant  gas  into  the  rib  region  and  improve  the 
overall  mass  transport.  This  demonstrates  the  importance  of 
accurate  permeability  values  used  in  the  numerical  simulation 
for  the  rib  (compressed)  portion  of  the  GDL. 

(4)  The  cross-channel  flow  also  affects  the  current  density  dis¬ 
tribution.  For  the  cases  tested,  it  is  found  that  the  peak 
current  density  occurs  in  the  high-pressure  channel  towards 
the  rib/channel  boundary  near  the  channel  with  a  low  pressure. 
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